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Abstract 

The set of bistochastic or doubly stochastic N x N matrices form a convex 
set called Birkhoff 's polytope, that we describe in some detail. Our problem 
is to characterize the set of unistochastic matrices as a subset of Birkhoff 's 
polytope. For iV = 3 we present fairly complete results. For = 4 partial 
results are obtained. An interesting difference between the two cases is that 
there is a ball of unistochastic matrices around the van der Waerden matrix 
for = 3, while this is not the case for = 4. 
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1. Introduction 

There is a surprising variety of contexts in which unistochastic matrices arise, 
and any one of them may be taken as a motivation for the present study. But 
let us first define our terms: An N x N matrix B is said to be bistochastic 
if its matrix elements obey 

i: Bij > u:J2 ^ij = 1 iii: J] A, = 1 ■ (1) 

i 3 

The name "bistochastic" has to do with the fact that these matrices are usu- 
ally supposed to act on probability distributions, thought of as component 
vectors. The first condition ensures that positive vectors arc transformed to 
positive vectors, the second that the sum of the components of the vectors re- 
mains invariant, and the third that the uniform distribution (a vector with all 
components equal) is a fixed point of the map. Hence a bistochastic matrix 
causes a kind of contraction of the probability simplex towards the uniform 
distribution. Condition iii is important because it guarantees that the map 
increases entropy. We obtain a bistochastic matrix if we start with a unitary 
matrix U and take the absolute value squared of its matrix elements, 

Bij = |L/jj|2 . (2) 

For connaisseurs of linear algebra, B is the Hadamard product of U and its 
complex conjugate. If there exists such a U then B is said to be unistochas- 
tic. If U is also real, that is orthogonal, then B is said to be orthostochastic. 
Bistochastic matrices arise frequently in situations where probability distri- 
butions are changing, and we will soon see why one may want them to be 
unistochastic. 

A somewhat distinguished bistochastic matrix is the van der Waerden 
matrix 5^, whose matrix elements obey 

The van der Waerden matrix is unistochastic. A corresponding unitary ma- 
trix is known as a complex Hadamard matrix. An example of a complex 
Hadamard matrix is the Fourier matrix, whose matrix elements are 
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U,, = ^q'^, 0<z,j<N-l. (4) 
V iv 



Here q = e^'^'/^ is a root of unity. Complex Hadamard matrices have a long 
history in mathematics [1]. If the matrix U is also real it is referred to simply 
as a Hadamard matrix. By the way. the name of the van der Waerden matrix 
has to do with the conjecture that this matrix has the largest permanent of 
all bistochastic matrices; this is true but took a long time to prove [2]. 

It is clear that we now have two mathematical questions on our hands: 

I: Given a bistochastic matrix, is it unistochastic? 

II: If so, to what extent is U determined by Bl 
The answer to question I turns out to depend on the bistochastic matrix 
chosen, so that in effect question I turns into the problem of characterizing 
the unistochastic subset of the set of all bistochastic matrices. But why are 
these questions interesting? The answer is that they naturally turn up in 
many contexts. Let us give a partial list of those. 

The first context has to do with the foundations of quantum mechanics. 
Here there are a number of approaches where one begins by arguing that 
transition probabilities, suitably defined, form bistochastic matrices. In at- 
tempting to build some group structure into these transition probabilities 
one is then led to require that they form unistochastic matrices, and so one 
runs into question I. A sample of the literature includes Lande [3], Rovelli 
[4] and Khrennikov [5]. 

The second context is classical computer science, especially the theory of 
error correcting codes, design theory, and other areas of discrete mathematics 
where real Hadamard matrices are very useful. Hadamard conjectured that 
such matrices exist when N = 2 and N = Ak [6]. The conjecture is still 
open, although much is known [7]. For explicit examples of real Hadamard 
matrices of all orders up to 256 x 256, consult Sloane [8]. 

The third context is quantum information theory, where the restriction 
to real Hadamard matrices is somewhat unnatural. Complex Hadamard 
matrices have been studied in the quantum optics community in the guise 
of symmetric multiports; they are examples of specially designed unitary 
transformations that can be reahzed in the laboratory [9] [10]. There is 
also an interesting theorem [11] to the effect that the classification of all 
possible teleportation schemes can be reduced to the classification of all sets 
of maximally entangled bases; this is relevant here because such sets can be 
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obtained from the combination of a Latin square and a complex Hadamard 
matrix. The construction of all possible Latin squares has nothing to do with 
us here, but the construction of all complex Hadamard matrices certainly has. 
Mathematicians have studied this problem with various motivations [12] [13] 

[14] [15]. 

The fourth context is the attempt to formulate quantum mechanics on 
graphs (in the laboratory on thin strips of, say, gold film). Here question 
I arises as a question about what Markov processes that have a quantum 
counterpart in the given setting [16] [17] [18]. In this connection studies of 
the spectra and entropies of unistochastic matrices chosen at random have 
been made [19]; we will return to some of these issues below. 

A fifth context is particle physics, where the interest centers on ques- 
tion II. Thus in the theory of weak interactions we encounter the unitary 
Kobayashi-Maskawa matrices (one for quarks and one for neutrinos), and 
Jarlskog raised the question to what extent such a matrix can be parametrized 
by the easily measured moduli of its matrix elements. The physically inter- 
esting case here is = 3 [20], and possibly also = 4 in case a fourth 
generation of quarks should be discovered [21]. (The same question occurs 
in scattering theory, with no restriction on N [22].) 

We end the list of possible applications here, and turn to the organisation 
of our paper. In section 2 we describe the set of all bistochastic matrices Bn- 
It is a convex polytope well known to practioneers of linear programming; it 
is sometimes called the assignment polytope because it arises in the problem 
of assigning A" workers to A^ tasks, given their efficiency ratings for each 
task. We describe the cases A^ = 3 and A^ = 4 in detail (A^ = 2 is trivial). In 
section 3 we discuss some generalities concerning unistochastic matrices, and 
then we characterize the unistochastic subset for the case of A^ = 3. Most of 
our results can be found elsewhere but, we believe, not in this coherent form. 
In section 4 we address the same question for A^ = 4. This turns out to be a 
more difficult task, but at the end of this section there will be a proof that 
every neighbourhood of the van der Waerden matrix contains matrices that 
are not unistochastic. This is a striking difference to the A^ = 3 case. Section 
5 summarises our conclusions. Some technical matters are found in three 
appendices. Our results on A" > 4 will be reported in a separate publication 
[23]. 



4 



2. Birkhoff's polytope 

The set Bn of bistochastic N x N matrices has {N — 1)^ dimensions. To 
do this count, note that the last row and the last column are fixed by the 
conditions that the row and column sums should equal one. The remaining 
(A^ — 1)^ entries can be chosen freely, within limits. Birkhoff proved that Bn 
is a convex polytope whose extreme points, or corners, arc the N\ permuta- 
tion matrices [24]. All corners are equivalent in the sense that they can be 
taken into each other by means of orthogonal transformations. A bistochas- 
tic matrix belongs to the boundary of B^ if and only if one of its entries is 
zero. The boundary consists of corners, edges, faces, 3-faccs and so on; the 
highest dimensional faces are called facets and consist of matrices with only 
one zero entry. For a detailed account of Bn, especially its face structure, see 
Brualdi et al. [25]. We will be even more detailed concerning B^ and ^4. For 
definiteness all 24 permutation matrices that occur when N — 4 are listed in 
Appendix A. 

It is convenient to regard the convex polytope Bn as a subset of a vector 
space with the van der Waerden matrix as its origin. The distance squared 
between two matrices is chosen to be 

D^{A,B) = Tr{A- B){A^ - B^) . (5) 

The distance squared between an arbitrary bistochastic matrix B and B^, is 
then given by 

= (6) 

111 particular, the distance between and a corner of the polytope becomes 
D = — 1. Permutations of rows or columns are orthogonal transforma- 
tions since they preserve distance and leave 5^ invariant. They also take 
permutation matrices (corners) into permutation matrices, hence they are 
symmetry operations of Birkhoff's polytope as well. 

The (Shannon) entropy of a bistochastic matrix is defined as the entropy 
of the rows averaged over the columns, 

i j 



5 



Its maximum value In N is attained at B^,. For some of its properties consult 
Slomczynski [26] et al. [19]. 

When N — 2 there are just two permutation matrices and B2 is a line 
segment between these two points. A general bistochastic matrix can be 
parametrized as 

c = cos^, s = sin^, 0<e<-. (8) 

When = 3 we have six permutation matrices forming the vertices of a four 
dimensional polytope. It admits a simple description: 




Theorem 1 : The 6 corners of are the corners of two equilateral triangles 
placed in two orthogonal 2-planes and centered at B^. 

The proof is easy, using as corners the permutation matrices Pq, Pi, . . . , P5 
from Appendix A. The two equilateral triangles are the convex combinations 



Al = poPo + psPa + PiPi 



PO P3 Pi 

Pa Po P3 
P3 Pa Po 



P0+P3+ Pa 



(9) 



and 



A2 = piPi + P2P2 + P5P5 



Pi P2 Pb 
P2 P5 Pi 
P5 Pi P2 



Pi + P2 + P5 = 1 • (10) 



The calculation we have to do is to check that D^^PqjP^) = D'^[Pq,P:i) — 
D^^Ps, P4) = 6 and similarly for the other triangle, and also that 



Tr(Ai-SO(A^-St) = 



(11) 



for all values of pi. This is so. 

There are thus 6 corners and 6-5/2 = 15 edges, all of which are extremal. 
The last is a rather exceptional property; in 3 dimensions only the simplex 
has it. There are 9 short edges of length squared = 4 and 6 long edges 
of length squared = 6, namely the sides of the two equilateral triangles. 
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Figure 1: Left: Birkhoff's polytopc for iV = 2 (centered at B^). Right: The 
graph of Birkhoff's polytope for = 3; single hnes have = 4 and double 
= 6. The double edges form the triangles mentioned in Theorem 1. 

A useful overview of is given by its graph, where we exhibit all corners 
and all edges (see fig. 1). All the 2- faces are triangles with one long and two 
short edges. The 3-faces in a 4 dimensional polytope are called facets and 
are made of matrices with a single zero. They are irregular tetrahedra with 
two long edges, one from each equilateral triangle (see fig. 4). 

The volume of is readily computed because it can be triangulated 
using only three simpfices. The total volume is 9/8. As N grows the total 
volume of Bn becomes increasingly hard to compute; mathematicians know 
it for iV < 10 [27]. 

The next case is the 9 dimensional polytopc B4. It has 24 corners and 276 
edges. The latter come in four types and we give the classification includ- 
ing the angle they subtend at and whether they consist of unistochastic 
matrices or not (see sections 3 and 4): 




Length squared Unistochastic Angle at origin Number of edges 

4C/ 4 Yes Acute 72 

6 6 No 90 degrees 96 

8 8 No Obtuse 72 

8U 8 Yes Obtuse 36 

All edges except the 8U ones are extremal. The 2-faces consist of triangles 
and squares. (Interestingly, for all it is true that the 2-faces of Birkhoff's 
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polytope Bn arc cither triangles or rectangles [25].) There are 18 squares 
bounded by edges of type 41/ and their diagonals are of type 8U. Three 
squares meet at each corner. If we pick four permutation matrices we obtain 
a 3-face, with six exceptions. The exceptions form 6 regular tetrahedra cen- 
tered at i?^ whose edges are non-extremal 8U edges. They are denoted Tj 
and explicitly listed in Appendix A; an example is 



Ti = po-Po + P7P7 + P16-P16 + P23-P23 = 



Po 


P7 


Pw 


P23 


P7 


Po 


P23 


Pw 


PI6 


P23 


Po 


P7 


P23 


PI6 


P7 


Po 



(12) 



When regular tetrahedra are mentioned below it is understood that we refer 
to one of these six. In a sense the structure can now be drawn; see fig. 2. 
The facets consist of matrices with one zero, so there are 16 facets. 

A subset of B4 that has no counterpart for B3 is the set of matrices that 
are tensor products of two by two bistochastic matrices. This subset splits 
naturally into several two dimensional components, and it turns out that 
they sit in B4 as doubly ruled surfaces inside the regular tetrahedra. Thus 
the following matrix, parametrised with two angles, is a tensor product of 
two matrices of the form (8): 



„2 2 
^1^2 


„2 2 
'"1*2 


2 2 
S1C2 


„2 2 
*1*2 


^2 2 
^1*2 


^2 2 


„2 2 
*1*2 


^2 2 
*l''2 


„2 2 


„2 2 


„2 2 


„2 2 




*1*2 


C1C2 


'"1*2 


„2 2 








*1*2 


^2^2 
*l''2 


^2 2 


^2 2 



ci = COS 6*1 etc. 



(13) 



These matrices form a doubly ruled surface inside the regular tetrahedron 
(12), analogous to that depicted in fig. 4. 

An interesting way to view B4, and one that will recur in section 4, stems 
from the following observation: 



Theorem 2: The 24 corners of B4 belong to a set of nine orthogonal hy- 
perplanes through Each regular tetrahedron belongs to six hyperplanes 
and contains the normal vectors of the remaining three hyperplanes. Each 
hyperplane contains four regular tetrahedra and its normal vector is the in- 
tersection of the remaining two regular tetrahedra. 
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Figure 2: How to begin to draw the surface of B4. Two tetrahedra whose 
edges are the non-extremal diagonals of squares are shown. The dashed line 
goes through the polytope; it connects the midpoints of two opposing 8U 
edges of two tetrahedra that are otherwise disjoint. 




Figure 3: A regular tetrahedron centered at B^. It contains the normal 
vectors of three orthogonal hypcrplanes and belongs entirely to another six. 
There are six such regular tetrahedra and pairs of them intersect along the 
normal vectors they contain. (Note that the dashed hne in Fig. 2 represents 
such a normal vector.) 
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Again the proof is a simple calculation, once the explicit form of the hy- 
perplanes is known. They are denoted Hj and listed in Appendix A. Prom 
now on, "hyperplane" always refers to one of these nine. Fig. 3 in a sense 
illustrates the theorem. 

It is quite helpful to have an incidence table for tetrahedra and hyper- 
planes available. It is 





Hi 


n2 


n3 


Hi 


n5 


He 






ng 






X 


X 


X 




X 


X 


X 




T2 




X 


X 


X 


X 




X 




X 




X 




X 




X 


X 


X 


X 






X 


X 






X 


X 


X 




X 


n 


X 




X 


X 


X 






X 


X 




X 


X 




X 




X 




X 


X 



(14) 



where the tetrahedra Tj and the hyperplanes Ilj are listed in Appendix A. 

For later purposes we will need some information about exactly how the 
hyperplanes divide the space into 2^ hyperoctants. For this reason we look 
at the rays 

Bi{t)^B^ + tVi , (15) 

where 1^ is a vector constructed in terms of the normal vectors rii, . . . , ng of 
the hyperplanes (see Appendix B), namely 
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-3 


-3 


-3 


1 


-3 


1 


1 


1 


4 


-3 


1 


1 


1 




-3 


1 


1 


1 



V2 = ni+n2+n3+n4+n5+ne+nr+ns—ng = - 



7 
-1 
-1 
-5 



-1 
-1 
-1 
3 



-1 
-1 
-1 
3 



-5 
3 
3 

-1 
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V3 = ni + n2 + n3 + n4 — n5 + n6 + n7 + ns — ng^ - 



5 
1 

-3 
-3 



1 - 
3 

1 - 
1 



3 
1 
3 
5 



3 
1 
5 
3 



( 
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All other cases can be obtained from one of these three by permutations of 
rows and columns. The various hyperoctants are convex cones centered on 
these rays. This gives a classification of the hyperoctants into six different 
types (since the parameter t can be positive or negative) called respectively 
type I±, II± and III±. Type I has 16 representatives and is especially note- 
worthy. For type I_ the centered ray hits the boundary in the center of one 
of the 16 facets, at the matrix |). In the other direction we also hit 

quite distinguished points. There are 16 ways of setting one entry of a bis- 
tochastic matrix equal to one, and this gives rise to 16 copies of B3 sitting in 
the boundary of B4. For the octants 1+ the centered ray hits the boundary 
precisely at the center of such a. B3, at the matrix 

In section 4 we will see how the structure of the unistochastic subset 
is related to the structure of Birkhoff 's polytope, and in particular to the 
features we have stressed. 
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3. The unistochastic subset, mostly N = 3 



Let us begin with some generahties concerning the unistochastic subset Un 
of Bn- The dimension of Bn is {N — 1)^ and the dimension of U{N) is 
A^^. Therefore the map U{N) Bn cannot be one-to-one. Now it is clear 
that muhiplying a row or a column by a phase factor — an operation that we 
refer to as rephasing — will result in the same bistochastic matrix via eq. (2). 
Therefore the map is naturally defined as a map from a double coset space 
to Bn- The double coset space is 

f/(l) X ■ ■ ■ X f/(l) \ U{N) I U{\) X ■ ■ ■ X f/(l) , (19) 

with iV f/(l) factors on the right and — 1 factors on the left, say. The 
dimension of this set is {N — 1)^ so now the dimensions match. There is a 
complication because the double coset space is not a smooth manifold. The 
action from the left of the C/(l) factors on the right coset space (in itself a well 
behaved flag manifold) has fixed points. These fixed points are easy to locate 
however (and always map to the boundary of Bn\ so that for most practical 
purposes we can think of our map as a map between smooth manifolds. 

In general we will see that the image of our map is a proper subset of ^Bjv, 
and the map is many-to-one. There is not much we can usefully say about 
the general case, except for two remarks: The unistochastic subset Un has 
the full dimension (A^ — 1)^ while the unistochastic subset of the boundary 
of Bn has dimension (A^— 1)^ — 2; why this is so will presently become clear. 

For N = 2 every bistochastic matrix is orthostochastic. A unitary matrix 
that maps to the matrix in eq. (8) is 



U 



c s 
s —c 



c = cose, s = sin^, 0<e<-. (20) 

2 



The matrix is given in dephased form. This means that the first row and the 
first column is real and positive. This fixes the U{1) factors mentioned above 
(unless there is a zero entry in one of these places) and from now on we will 
present all unitary matrices in this form. For any A^ it is straightforward 
to check whether a given edge of Bn is unistochastic. For A^ = 3 the edges 
of length squared equal to 4 are unistochastic, and for A^ = 4 we have the 
results given in table (12). 
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Given a 3 x 3 bistochastic matrix it is easy to check whether it is unis 
tochastic or not [28 
the matrix 



We form the moduh Vij 



Bij and write down 



U = 



^20 



^216 



8021 



(21) 



If this matrix is unitary the original matrix is unistochastic. The unitarity 
conditions simply say that the first two columns are orthogonal; the last 
cohimn by construction has the right moduli and does not impose any further 
restrictions. Therefore the problem is to form a triangle from three line 
segments of given lengths 

Lo = rooToi Li = riorii L2 = r2or2i . (22) 
This is possible if and only if the "chain-hnks" conditions are fulfilled, i.e. 



\Li — L2I < Lq < ivi + L2 



(23) 



The bistochastic matrix corresponding to U sits at the boundary of U3 if 
and only if one of these inequahties is saturated. When eq. (23) holds the 
solution is 



cos 011 



LI 



t2 _ t2 
^0 ^1 



cos (011 - 02l) 



COS 021 



^2 



2i^0-^2 



-^0 -^1 ^2 



(24) 



(25) 



2L1L2 

There is a two-fold ambiguity (corresponding to taking the complex conjugate 
of the matrix, U — > U*). The area A of the triangle is easily computed and 
the chain-links conditions are equivalent to the single inequality A > 0. As 
a matter of fact we can form six so called unitarity triangles in this way, 
depending on what pair of columns or rows that we choose. Although their 
shapes differ their area is the same, by unitarity [20] . 

Because we can easily decide if a given matrix is unistochastic, it is easy 
to characterize the unistochastic set U3. We single out the following facts 
(some of which are known [28]) for attention: 



13 



Figure 4: Birkhoff's polytope for = 3. Left: One of the two orthogo- 
nal equilateral triangles centered at B^, with its Tinistochastic subset (the 
boundary is the famous hypocycloid). Right: A facet, an irregular tetrahe- 
dron, with its doubly ruled surface of unistochastic matrices. 

Theorem 3 : The unistochastic subset U3 of B3 is a non-convex star shaped 
four dimensional set whose boundary consists of the set of orthostochastic 
matrices. It contains a unistochastic ball of maximal radius \/2/3, centered 
at Bi,. The set meets the boundary of B3 in a doubly ruled surface in each 
facet. 

The relative volume of the unistochastic subset is, according to our numerics, 

^ 0.7520 ± 0.0005 . (26) 

We did not attempt an analytical calculation; details of our numerics are in 
Appendix B. 

Theorem 3 is easy to prove. To see that U3 is non-convex we just draw its 
intersection with one of the equilateral triangles that went into the definition 
of the polytope, and look at it (see fig. 4). An amusing side remark is 
that the boundary of the unistochastic set in this picture is a 3-hypocycloid 
[19]. It can be obtained by roUing a circle of radius 1/3 inside the unit 
circle. The maximal unistochastic ball is centered at B^, and touches the 
boundary at the hypocycloid, as one might guess from the picture; its radius 
was deduced from results presented in ref. [29]. To see that the boundary 
consists of orthostochastic matrices is the observation that when the chain- 
links conditions are saturated the phases in U will equal ±1. That the 



14 



set is star shaped then follows from an explicit check that there is only 
one orthostochastic matrix on any ray from 5^. Finally fig. 4 includes an 
exphcit picture of the unistochastic subset of a facet. The reason why it has 
codimension one is that a matrix on the boundary of Bn has a zero entry, 
which means that the number of phases available in the dephased unitary 
matrix drops with one, and then the dimension of the unistochastic set also 
drops with one; the argument goes through for any A^. 

Finally let us make some remarks on entropy. We compare the Shannon 
entropy averaged over B3 using the flat measure, the Shannon entropy aver- 
aged over W3 also using the flat measure, and the maximal Shannon entropy 
'S'max- Numerically we find that 



with all digits significant. Observe that the latter average is larger since some 
matrices of small entropy close to the boundary of B3 are not unistochastic 
and do not contribute to the average over U3. The above data may be 
compared with the maximal possible entropy S'max = ln3 ~ 1.099, attained 
at Bi,, and also with 



which is the average taken over W3 with respect to measure induced by the 
Haar measure on U{3). This analytical result follows directly from the work of 
Jones, who computed the average entropy of squared components of complex 
random vectors [30]. It is easy to see that the two averages coincide. For 
details of our numerics consult Appendix B. 



(S)b3 ^ 0.883 and {S)us ^ 0.908 



(27) 




1 + - « 0.833 , 

2 3 



(28) 



15 



4. The unistochastic subset, mostly N = 4 



The case N — A is more difficult. It is also clear from the outset that it 
will be qualitatively different — thus the dimension of the orthogonal group 
is too small for the boundary of the unistochastic set W4 to be formed by 
orthostochastic matrices alone. There are other differences too, as we will 
see. 

Given a bistochastic matrix we can again define rij — and consider 



roo 




ro2 


• 


rio 






• 


r2o 






• 








• 



(29) 



Phases must now be chosen so that this matrix is unitary, and more especially 
so that the three columns we focus on are orthogonal. Geometrically this is 
the problem of forming three quadrilaterals with their sides given and six 
free angles. This is not a simple problem, and in practice we have to resort 
to numerics to see whether a given bistochastic matrix is unistochastic (see 
Appendix B for details). There are some easy special cases though. One 
easy case is that of a matrix belonging to the boundary of Bn- Then the 
matrix U must contain one zero entry and when we check the orthogonality 
of our three columns two of the equations reduce to the problem of forming 
triangles, so that the angles are completely fixed when we consider the final 
orthogonality relation. Another easy case concerns the regular tetrahedra. 
They turn out to consist of orthostochastic matrices; for the example given 
in eq. (12) a corresponding orthonormal matrix is 



(30) 



/PO VP7 sJV\% VP23 

VPw VP23 -y/Po 
VP23 -y/Pl6 -VPO . 

This saturates a bound saying that the maximum number oi N x N permu- 
tation matrices whose convex hull is unistochastic is not larger than 2!^], 
where [A^/2] denotes the integer part of N/2 [31]. 

Let us now turn our attention to B^. Hadamard [6] observed that up 
to permutations of rows and columns the most general form of the complex 
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Hadamard matrix is 



1111 

1 e**^ —1 —e^'^ 
1-11-1 
1 -e'"^ -1 e"^ 



(31) 



One can show that this is a geodesic in U{N). The news, compared to 
= 3, is that is orthostochastic because H{0) is reaL Moreover there 
is a continuous set of dephased unitaries mapping to the same B. In a 
calculational tour de force, Auberson et al. [21] were able to determine all 
bistochastic matrices whose dephased unitary preimages contain a continuous 
ambiguity (and they found that the ambiguity is given by one parameter in 
all cases). There are three such famihes. Using the notation of ref. [21] they 
consist of matrices of the following form: 



Type A: 



abed 

bade 
e f 9 h 
f e h g 



Type B: 



Type C: 



a a 

b b 

e e 

d d 



a 
b 

e 
d 



l-a 
\-d 



(32) 



„2 2 
*1*2 


^2„2 


^2„2 
•--3 "^2 


^2 2 


„2 2 
^1^2 


„2 2 
C1C2 


„2 2 

< ;v^2 


„2 2 
"'3"'2 


„2 2 


„2 2 


„2 2 


„2 2 
C354 


„2 2 


„2 2 


„2 2 
63C4 


„2 2 
C3C4 



(33) 



Here Ci = cos^^i, Si = sin^i, and so on. Type A consists of nine five di- 
mensional sets, type B of nine four dimensional sets, and type C of six three 
dimensional sets. In trying to understand their location in B4 the observa- 
tion in section 2 concerning the nine orthogonal hyperplanes begins to pay 
dividends. (In particular, consult the incidence table 14.) Type A consists of 
the linear subspaces obtained by taking all intersections of four hyperplanes 
that contain exactly two regular tetrahedra. Type C consists of the linear 
subspaces obtained by taking all intersections of six hyperplanes that contain 
no permutation matrices at all. Type B finally consists of curved manifolds 
confined to one hyperplane. Auberson's families are not exclusive. In par- 
ticular tensor product matrices belong to families A and B, which means 
that there are two genuinely different ways of introducing a free phase in the 
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corresponding unitary matrix. Outside the three sets A, B and C Auberson 
et al. find a 12-fold discrete ambiguity in the dephased unitaries, dropping 
to 4-fold for symmetric matrices [21]. 

Tensor product matrices S4 = i?2 ® B'2 appear because 4 = 2 x 2 is a 
composite number. That they are always unistochastic follows from a more 
general result: 

Lemma 1 : Let Bk and Bm be unistochastic matrices of size K and M, 
respectively. Then the matrix Bj^ — B^ (8) Bm of size KM is unistochastic. 
The corresponding unitary matrices contain at least {K — 1)(M — 1) free 
phases when dephased. 

That Bjsr is unistochastic follows from properties of the Hadamard and the 
tensor products. By definition, the Hadamard product Ao B of two matrices 
is the matrix whose matrix elements arc the prodiicts of the corresponding 
matrix elements of A and B. Then Bk = Uk ° and Bm = Um ° Um 
implies that Bn = {Uk o U*k) ® {Um o U^^) = {Uk O Um) o {U*k ® U*m), so 
it is unistochastic. The existence of free phases is an easy generalization of 
proposition 2.9 in Haagerup [13]. 

The hyperplane structure of B4 reverberates in the structure of the unis- 
tochastic set in several ways. Let us consider how the tangent space of U {N) 
behaves under the map to Bn- In equations, this means that we fix a unitary 
matrix Uq and expand 

U{t) = e^'**C/o = (1 + iht - \hH'' + . . . )C/o (34) 

where h is an Hermitian matrix. Then we study bistochastic matrices with 
elements -By(i) = to first order in t. The following general features 

are observed: 

• Generically the tangent space of U {N) maps onto the tangent space of 
Bjv. This implies that the dimension of the unistochastic set is equal to 
that of Bn; we checked this statement by generating unitary matrices 
at random using the Haar measure on the group. 

• A matrix element in B receives a first order contribution only if it is 
non- vanishing. Hence the map of the tangent space of U{N) to the 
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tangent space of Bn is degenerate at the boundary of the polytope. 
In general such behaviour is to be expected at the boundary of the 
unistochastic set Un- 

• If Uq is real the map is degenerate in the sense that the tangent space 
maps to an A^(A^ — l)/2 dimensional subspace of the tangent space of 
Bn. 

• If Uq maps to a corner of the polytope then the first order contributions 
vanish. To second order we pick up the tip of a convex cone whose 
extreme rays are the N{N — 1)/2 4C/ edges emanating from that corner. 

For = 4 the story becomes interesting when we choose Uq equal to the 
Hadamard matrix H{(j)). Then we find that the tangent space at Uq maps 
into one of the nine hyperplanes; which particular one depends on how we 
permute rows and columns in eq. (31). The question therefore arises whether 
the orthostochastic van der Waerden matrix belongs to the boundary of the 
unistochastic set — or not since a priori such degeneracies can occur also in 
the interior of the set. 

We know that we can form curves of unistochastic matrices starting from 
and moving out into the nine hyperplanes. Can we form such curves that 
go directly out into one of the 2^ hyperoctants? Here the division of the 2^ 
hyperoctants into six different types becomes relevant. We have investigated 
whether their central rays given in eqs. (15-18) consist of unistochastic ma- 
trices, or not. Let us begin with the 16 hyperoctants of type I, where the 
central ray Bi(t) — + tVi hits the boundary in the center of one of the 
16 Bs sitting in the boundary (at t = 1/3), and in the center of one of the 
16 facets (at t = —1/9). Of these two points, the first is unistochastic, the 
second is not. A one parameter family of candidate unitary matrices that 
maps to the central ray is 



(35) 



where i > and we permuted the columns relative to eq. (16) in order to 
get the unitarity equations in a pleasant form. (We do not need to give the 
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phases for the last column.) The conditions that the first three columns be 
orthogonal read 



gi0n ^ gi<^2i ^ gi<^3i + L = (36) 

gi012 + gi<^22 + gi<^32 + X = (37) 
gi(011-012) _j_ gi(021-022) _j_ giC'/'si-te) -)- X/ = (38) 



where 



l-3t 



(39) 



1 + t 

In Appendix C we prove that the system of equations (36-38) 

1. has no real solutions for L > 1, 

2. for < L < 1 has the solution 

^^"1 > cos0=^ = -:^. (40) 

It follows that the central ray is unistochastic for the hyperoctants of type 

1+ (and the unitary matrices on the central ray tend to the real Hadamard 
matrix at t = 0). In the other direction the central ray is not unistochastic 
for type I_. Thus we have proved 



Theorem 4 : For N — A there are non-unistochastic matrices in every neigh- 
bourhood of the van der Waerden matrix 5^. At the map U{4) B4 
aligns the tangent space of U (4) with one of the nine orthogonal hyperplanes. 

The structure of the unistochastic set is dramatically different depending on 
whether = 3 or = 4. It is only in the former case that there is a ball of 

unistochastic matrices surrounding S^. On the other hand, the hyperoctants 
are not empty — some of them do contain unistochastic matrices all the way 
down to Bi,. 

Concerning the other hyperoctants, for types II_, III+, and III_ the cen- 
tral rays hit the boundary of the polytope in points that are not unistochastic, 
but numerically we find that a part of the ray close to is unistochastic. 
For type 11+ we hit the boundary in a unistochastic point and numerically 
we find the entire ray to be unistochastic. There is still much that we do 
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not know. We do not know if the hyperoctants of type I_ are entirely free 
of unistochastic matrices, nor do we know if W4 is star shaped, or what its 
relative volume may be. What is clear from the results that we do have is 
that the global structure of Birkhoff's polytope reverberates in the struc- 
ture of the unistochastic subset in an interesting way — it is a little bit like a 
nine dimensional snowflakc, because the nine hyperplanes in 84^ can be found 
through an analysis of the behaviour of ZY4 in the neighbourhood of B^. 
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5. Conclusions 

Our reasons for studying the unistochastic subset of Birkhoff 's polytope have 
been summarized in the introduction. Because the problem is a difficult 
one we concentrated on the cases = 3 and N = 4. Our descriptions of 
Birkhoff 's polytope for these two cases are given in Theorems 1 and 2, respec- 
tively, and a characterization sufficient for our purposes of the unistochastic 
set for = 3 was given in Theorem 3. For N — A the dimension of the unis- 
tochastic set is again equal to that of the polytope itself, but its structure 
differs dramatically from the = 3 case. In particular Theorem 4 states 
that for = 4 there arc non-unistochastic matrices in every neighbourhood 
of the van der Waerden matrix. Hence there does not exist a unistochastic 
ball surrounding the van der Waerden matrix. We observed that the struc- 
ture of the unistochastic set at the center of the polytope reflects the global 
structure of the latter in an interesting way. 

It is natural to ask to what extent the difference between the two cases 
reflects the fact that 3 is prime while 4 is not. Although this is not the place 
to discuss the cases A^ > 4, since some of us intend to do so in a separate 
publication [23], let us mention that the dimension of the unistochastic set 
is equal to that of for all values of A^. On the other hand it is only in the 
case of A^ being a prime number that we have been able to show that there 
is a unistochastic ball surrounding the van der Waerden matrix. 
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Appendix A: Notation 

An explicit list of permutation matrices for A'^ = 4 is 
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10 
10 
10 
1 



Pi 



10 
10 
1 
10 



10 
10 
10 
1 



(41) 



10 
10 
1 
10 



P4 = 



10 
1 
10 
10 



10 
1 
10 
10 



(42) 



10 
10 
10 
1 



10 
10 
1 
10 



10 
10 
10 
1 



(43) 



10 
10 
1 
10 



10 
1 
10 
10 



10 
1 
10 
10 

(44) 



10 

10 
10 
1 



13 



10 

10 
1 
10 



14 



10 

10 
10 
1 



(45) 



10 
10 

1 
10 



Pi 



16 



10 
1 

10 
10 



17 



10 
1 
10 
10 

(46) 



23 



18 



21 



1 

10 
10 
10 



1 
10 
10 
10 



Pl9 = 



22 



1 
10 
10 
10 



1 
10 
10 
10 



P20 — 



1 
10 
10 
10 



(47) 
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1 
10 
10 
10 



(48) 

In Birkhoff's polytope these 24 matrices form the corners of 6 regular tetra- 
hedra, namely the convex hulls of the sets 

— {Po, P7, Pl6, P23} T2 — {Pi, Pq, Pi7, P22} 

Ts = {P2, Pio: As, P21} T4 = {Pa, Pn, P12, P20} (49) 

^5 = {P4! Ps) Pl5! P19} Tq = {P5, Pg, Pi4, Pis} . 

The nine hyperplanes mentioned in Theorem 2 consist of matrices of the form 



n. 



Poo 


Poi • 


Pio 


Pll • 


• 


• • 


• 


• • 


Poo 


Poi • 


• 


• • 


• 


• • 


P30 


P31 • 


Poo 


• P02 


• 


• • 


P20 


• P22 


• 


• • 



Hp 



Poo 


Poi • 


• 


• • 


P20 


P21 • 


• 


• • 


Poo 


• P02 


Pio 


• P12 


• 


• • 


• 


• • 


Poo 


• P02 


• 


• • 


• 


• • 


P30 


• P32 



(50) 



(51) 



(52) 
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-Boo 


• 


• 


-Bo3 


-Bio 


• 


• 


-Bia 


• 


• 


• 


• 


• 


• 


• 


• 



He 



-Boo 


• 


• 


-Boa 


• 


• 


• 


• 


-B20 


• 


• 


-B23 


• 


• 


• 


• 



-Boo 


• 


• 


-B03 


• 


• 


• 


• 


• 


• 


• 


• 


-B30 


• 


• 


-B33 



(53) 



(54) 



where the matrix elements that are exphcitly written are assumed to sum to 
one (and similarly for the remaining three blocks taken separately). 
The normal vectors of these hyperplanes are the matrices 



1 


1 


-1 


-1 


1 


1 


-1 


-1 


-1 


-1 


1 


1 


-1 


-1 


1 


1 



and so on. 



(55) 



Appendix B: Numerics 

I. Average entropy in B3. 

To generate a random bistochastic matrix according to the flat measure 
on Bs C M'', we have drawn at random a point {x, y, z, t) in the 4-dimensional 
hypercube. It determines a minor of a A'^ = 3 matrix, and the remaining five 
elements of may be determined by the unit sum conditions in eq. (1). 
Condition i is fulfilled if the sums in both rows and both columns of the 
minor does not exceed unity, and the sum of all four elements is not smaller 
than one. If this was the case, the random matrix B3 was accepted to the 
ensemble of random bistochastic matrices. If additionally, the chain links 
condition (23) were satisfied, the matrix was accepted to the ensemble of 
unistochastic matrices, generated with respect to the flat measure on U3. 
The mean entropies, (27), were computed by taking an average over both 
ensembles consisting of 10'' random matrices, respectively. 
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II Numerical verification, whether a given bistochastic matrix B is unis- 
tochastic. 

We have performed a random walk in the space of unitary matrices. Start- 
ing from an arbitrary random initial point Uq we computed Bq = Uqo f/* and 
its distance to the analyzed matrix, Dq = D{Bq,B), as defined in (5). We 
were fixing a small parameter a ^ 0.1, generated a random Hermitian matrix 
H distributed according to the Gaussian unitary ensemble [32], and found a 
unitary perturbation V — exp{—aH). The matrix Un+i — VUn was accepted 
as a next point of the random trajectory, if the distance Dn+i was smaller 
than the previous one, D„. If a certain number (say 100) of random matrices 
V did not allow us to decrease the distance, we were reducing the angle a 
by half, to start a finer search. A single run was stopped if the distance D 
was smaller then e = 10~^ (numerical solution found), or a got smaller then 
a fixed cut off value (say amin = 10~^). In the latter case, the entire pro- 
cedure was repeated a hundred times, starting from various unitary random 
matrices Uo, generated according to the Haar measure on f/(4) [33]. The 
smallest distance -Dmin and the closest unistochastic matrix Bmin = Un°Un 
was recorded. 

To check the accuracy of the algorithm we constructed several random 
unistochastic matrices, B = U oU^ and verified that random walk procedure 
was giving their approximations with L>min < e. 

Appendix C: A system of equations 

In order to curtail a plethora of indices in Eqs. (36-38) and ease the 
subsequent notation let us introduce shorthands: ipj = (pji, ipj = —(pj2, 
j — 1,2,3. With that the system rewrites as 




(56) 
(57) 
(58) 



We shall prove the following: 



Lemma : The system of equations (56-58) 
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1. has no real solutions for L > 1, 

2. for < L < 1 has the solution 



V'l = -0 -02 = V^a = ' ^ t + 1 2 ^ ^ 

unique up to obvious permutations, 
3. for L = 0, 1 has continuous families of solutions. 

Indeed, each of the unimodal numbers e^'^'', A; = 1, 2, 3 is a root of: 

P(A) = (A-e*'^i)(A-e*'^^)(A-e*^«) (60) 
— A"^ — (e*'^^ + e*'''^ + e*'^^)A^ + (e**^'^^"'"'^^) + e'*^v'i+</'3) _|_ g«(¥'2+¥'3)^^ _ g(*'pi+¥'2+¥'3) 

= A^ + A^L - ALe'* - e'* = A2(A + L) - (1 + AL)e**, 

where $ = y^i + + v^a, and we used (56) and the reality of L. Thus each 
A = e*^^ {k = 1,2,3), fulfils: 

A'(A + L) = (l + AL)e'*. (61) 

Analogously, = e^'^^ (fc = 1, 2, 3), fulfils 

//'(// + L) = (l + /xL)e'*, (62) 

with = ■^i + + ■03- 

Observe now, that if A = e^'^'' and // = e*^* are solutions of (56-58) with 
the same number k, {k — 1, 2, 3) then, upon the same reasoning applied to 
(58), Xfx fulfils 

AV'(AyU + L) = (1 + A/xL)e'(*+*) . (63) 

Multiplying (61) by (62) and finally by (63) after exchanging its sides, we 
obtain, after division by A^/x^e*'^*"''"^^ ^ 0, 

(L + A)(L + /i)(LA/i + 1) = (LA + 1)(L// + 1)(L + A/x), (64) 

which, upon substitution A = e**^*, /i — e^'^'^ and putting everything on one 
side factorizes to 

L{L - l)(e*^'= - l)(e*'^'= - l)(e*('^'=+^'=) - 1) = 0, (65) 



27 



(any computer symbolic manipulation program can be helpful in revealing 
(65) from (64)). 

Hence, if L 7^ 0,1, then for each pair {(pk,ipk), k — 1,2,3, either: a) 
one of the angles is zero or b) they are opposite. The latter case can not 
occur for all three pairs since then e^^^i+'^i) + e*('^2+V2) + e'^^3+4>3) ^ ^ -L, 
hence at least one of (^^ or ipk equals zero. Up to unimportant permutations 
we can assume 933 = 0, but then, since e"^^ + 6**^^ + e*^^ = —L G M, we 
immediately get (pi — — </?2- This determines also all other angles (also up to 
some unimportant permutation) and wc end up with the solution announced 
in point 2. above as the only possibility, but such a solution exists only if 
L < 1. 

To prove 3. observe that 

1. for L = 0, 

(/?! = (/?, (p2^(p + 27r/3, (p3^(fi + 47r/3, (66) 
= ^2^^ + 2n/3, V3 = V' + 47r/3, (67) 

is a legitimate solution of (56-58) for arbitrary (p and ip, 

2. for L = 1 

(Pl=(p, (P2 = 7T, ip3 = ip + 7r (68) 
Ipi^ -ip + TT, 1p2 = TT, = (69) 

is a solution for an arbitrary ip. 
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